Risk based evaluations of competing agronomic climate adaptation strategies: Workflow example for rice planting strategies from crop growth simulation evidence
This notebook presents a novel approach for developing and targeting robust recommendations from crop simulation, long term experiments or multi-season survey data based on a risk optimization model.
For the computational model, Terry Hurley at University of Minnesota wrote the first Matlab code from which we have adapted to Octave and R.
To reproduce the analyses, one requires Octave installation. You can download Octave from here: https://octave.org/download.
2 Stylized example
In order to run this step, one needs to install Octave. The folder should contain the different functions that used for the optimization including CreateWTPBoundsbyCell and CreateTableData.
library(rio)ID=rep(1,1000)g=as.data.frame(g*1000)q=as.data.frame(q*1000)f=as.data.frame(f*1000)Rep=1:nrow(g)Area=rep(1,1000)QvsG=data.frame(ID,Rep,Area,q,g)export(QvsG,colNames=F,"QvsG.xlsx")QvsG=as.matrix(QvsG)colnames(QvsG) <-c(1,2,3,4,5)library(rhdf5)h5save(QvsG,file="QvsG.hdf")QvsF=data.frame(ID,Rep,Area,q,f)export(QvsF,colNames=F,"QvsF.xlsx")QvsF=as.matrix(QvsF)colnames(QvsF) <-c(1,2,3,4,5)h5save(QvsF,file="QvsF.hdf")FvsG=data.frame(ID,Rep,Area,f,g)export(FvsG,colNames=F,"FvsG.xlsx")FvsG=as.matrix(FvsG)colnames(FvsG) <-c(1,2,3,4,5)h5save(FvsG,file="FvsG.hdf")# Experiments -------------# G vs Q: Q as baseline -----------------GvsQ=data.frame(ID,Rep,Area,g,q)export(GvsQ,colNames=F,"GvsQ.xlsx")GvsQ=as.matrix(GvsQ)colnames(GvsQ) <-c(1,2,3,4,5)h5save(GvsQ,file="GvsQ.hdf")## F vs Q: Q as baseline ---------------FvsQ=data.frame(ID,Rep,Area,f,q)export(FvsQ,colNames=F,"FvsQ.xlsx")FvsQ=as.matrix(FvsQ)colnames(FvsQ) <-c(1,2,3,4,5)h5save(FvsQ,file="FvsQ.hdf")## G vs F: F as baseline ----------------GvsF=data.frame(ID,Rep,Area,g,f)export(GvsF,colNames=F,"GvsF.xlsx")GvsF=as.matrix(GvsF)colnames(GvsF) <-c(1,2,3,4,5)h5save(GvsF,file="GvsF.hdf")
2.3 Q vs G
The computational model consists of three main functions.
2.3.0.1 SOSDIntegralTest function
function [result] = SOSDIntegralTestv3(DOMINANT,DOMINATED)
% Tests for second-order stochastic dominance assuming uniform distribution
% over vector elements
% Algorithm Adapted from Levy, H. (2006). Stochastic Dominance: Investment Decision
% Making Under Uncertainty (Second Edition). Springer. New York, NY.
% Pages 180-182.
% result = 1 implies suspected DOMINANT distribution is in fact SOSD over DOMINATED
% result = 0 implies suspected DOMINANT distribution is not SOSD over DOMINATED
len1 = length(DOMINANT);
len2 = length(DOMINATED);
if len1 ~= len2
error('Lengths do not match!');
end
DOMINANTSORT = sortrows(DOMINANT,1);
DOMINATEDSORT = sortrows(DOMINATED,1);
work = zeros(len1,2);
flag1 = 1; % stays 1 if cumulative of suspected dominated >= cumulative of suspected dominant
flag2 = 0; % stays 0 unless cumulative of supected dominated > cumulative suspected dominant for some value
for ind = 1:len1 % Check to see if suspected dominant distribution satisfies SOSD integral condition
% relative to the suspected dominated distribution
if ind == 1
work(ind,1) = DOMINANTSORT(ind,1);
work(ind,2) = DOMINATEDSORT(ind,1);
else
work(ind,1) = DOMINANTSORT(ind,1) + work(ind - 1,1);
work(ind,2) = DOMINATEDSORT(ind,1) + work(ind - 1,2);
end
if work(ind,1) < work(ind,2)
flag1 = 0; % Not SOSD because suspected dominant distribution has higher cumulative value than suspected dominated
end
if work(ind,1) > work(ind,2)
flag2 = 1; % To be SOSD the suspected dominated distribution must have higher cumulative
% distribution than suspected dominant distribution for at least one value
end
end
result = flag1 * flag2;
2.3.0.2 ConstRiskThresh function
function [Result] = ConstRiskThresh(Data, CompVarNum, BaseVarNum)
% INPUT Matrices: ndataseq01060216 and ndataseq10060216 from crop
% simulation models with Columns
% 1 = cell30m/Unique Cell ID
% 2 = Replication
% 3 = Maize Area in Cell (ha)
% 4 = Maize Yield (kg/ha) for CM1510 with 40 kg N / yldCM1040
% 5 = Maize Yield (kg/ha) for CM1510 with 0 kg N / yldCM1000
CompManNum = 1;
if CompVarNum == 5
CompManNum = 0;
end
BaseManNum = 1;
if BaseVarNum == 4
BaseManNum = 0;
end
% Outputs
% 1 = cell30m
% # = Comparison Cultivar ID (CM1510 = 10, CM1509 = 9, Improved = 0)
% 2 = Comparison Nitrogen ID
% # = Base Cultivar ID (CM1510 = 10, CM1509 = 9, Improved = 0)
% 3 = Base Nitrogen Management ID
% 4 = Mean Yield for Comp
% 5 = Standard Deviation of Yield for Comp
% 6 = CV of Yield for Comp
% 7 = Maximum Yield for Comp
% 8 = Minimum Yield for Comp
% 9 = Probability of Crop Failure for Comp
% 10 = Min Proportion for Comp to SOSD Base
% 11 = Mean Yield for Base
% 12 = Standard Deviation of Yield for Base
% 13 = CV of Yield for Base
% 14 = Maximum Yield for Base
% 15 = Minimum Yield for Base
% 16 = Probability of Crop Failure for Base
% 17 = Min Proportion for Base to SOSD Comp
% 18 = Difference in mean Comp - Base
% 19 = Difference in standard deviation Comp - Base
% 20 = Difference in CV Comp - Base
% 21 = Difference in Prob of Crop Failure Comp - Base
% 22 = Min Proportion for Comp to SOSD Base divided by average base yield
% 23 = Min Proportion for Base to SOSD Comp divided by average base yield
% 24 = Comp More Risky (-1)/less Risky (1)/ Indeterminant (0) compared to Base
% 25 = Maize Area
% Get list of cell30m
CELLIDS = Data;
CELLIDScondy = CELLIDS(:,2) ~= 1;
CELLIDS(CELLIDScondy,:) = [];
CELLIDS(:,2:3) = [];
LEN = length(CELLIDS);
Result = ones(LEN, 25);
count = 1;
%while count <= LEN
id = CELLIDS(count,1);
Yields = Data;
cellcond = Yields(:,1) ~= id;
Yields(cellcond,:) = [];
TEMP = SOSDConstBoundsv3(Yields(:,CompVarNum),Yields(:,BaseVarNum));
Result(count, 1) = id;
% Result(count, 2) = CompVar;
Result(count, 2) = CompManNum;
% Result(count, 4) = BaseVar;
Result(count, 3) = BaseManNum;
Result(count, 4:24) = TEMP.';
Result(count, 25) = Yields(1,3);
% count = count + 1;
%end
2.3.0.3 SOSDConstBounds
function [Result] = SOSDConstBoundsv3(DIST1,DIST2)
% Finds the min number where DIST1 (the comparison distribution) SOSD DIST2 (the base distribution) and
% the min number where DIST2 SOSD DIST1 in addition to a variety of
% descriptive statistics for DIST1 and DIST2
if size(DIST1,2) > 1
error('DIST1 must be a column vector!')
end
if size(DIST2,2) > 1
error('DIST2 must be a column vector!')
end
% Output Result row definitions
MEANDIST1 = 1;
SDDIST1 = 2;
CVDIST1 = 3;
MAXDIST1 = 4;
MINDIST1 = 5;
PRCROPFAILDIST1 = 6;
DIST1MINPROP = 7;
MEANDIST2 = 8;
SDDIST2 = 9;
CVDIST2 = 10;
MAXDIST2 = 11;
MINDIST2 = 12;
PRCROPFAILDIST2 = 13;
DIST2MINPROP = 14;
MEANDIFF = 15;
SDDIFF = 16;
CVDIFF = 17;
CFDIFF = 18;
RELCOMPSOSDBASE = 19;
RELBASESOSDCOMP = 20;
DELTARISK = 21;
Result = -999999 * ones(21,1); % Initialize Results
% Calculate the mean, standard deviation, coeffcient of variation, minimum
% and maximum for DIST1
Result(MEANDIST1,1) = mean(DIST1);
Result(SDDIST1,1) = std(DIST1);
if Result(MEANDIST1,1) > 0
Result(CVDIST1,1) = Result(SDDIST1,1) / Result(MEANDIST1,1);
end
Result(MAXDIST1,1) = max(DIST1);
Result(MINDIST1,1) = min(DIST1);
ldist1 = length(DIST1);
% Calculate the mean, standard deviation, coeffcient of variation, minimum
% and maximum for DIST2
Result(MEANDIST2,1) = mean(DIST2);
Result(SDDIST2,1) = std(DIST2);
if Result(MEANDIST2,1) > 0
Result(CVDIST2,1) = Result(SDDIST2,1) / Result(MEANDIST2,1);
end
Result(MAXDIST2,1) = max(DIST2);
Result(MINDIST2,1) = min(DIST2);
ldist2 = length(DIST2);
% Calculate differences in descriptive statistics between DIST2 and DIST1
if Result(MEANDIST1,1) == 0 && Result(MEANDIST2,1) == 0
Result(MEANDIFF,1) = -999999;
else
Result(MEANDIFF,1) = Result(MEANDIST1,1) - Result(MEANDIST2,1);
end
if Result(SDDIST1,1) == 0 && Result(SDDIST2,1) == 0
Result(SDDIFF,1) = -999999;
else
Result(SDDIFF,1) = Result(SDDIST1,1) - Result(SDDIST2,1);
end
if Result(CVDIST1,1) >= 0 && Result(CVDIST2,1) >= 0
Result(CVDIFF,1) = Result(CVDIST1,1) - Result(CVDIST2,1);
else
Result(CVDIFF,1) = -999999;
end
% Initialize Thresholds for Stopping Golden Section Search
thresh = 0.00001;
bailthresh = 1000000;
if Result(MEANDIST1,1) > 0 || Result(MEANDIST2,1) > 0
% Calculate Probability of crop failiure for DIST1
CropFailDIST1 = DIST1;
CropFailDIST1cond = DIST1(:,1) ~= 0;
CropFailDIST1(CropFailDIST1cond,:) = [];
Result(PRCROPFAILDIST1,1) = length(CropFailDIST1) / ldist1;
% Calculate Probability of crop failiure for DIST1
CropFailDIST2 = DIST2;
CropFailDIST2cond = DIST2(:,1) ~= 0;
CropFailDIST2(CropFailDIST2cond,:) = [];
Result(PRCROPFAILDIST2,1) = length(CropFailDIST2) / ldist2;
tcomp = Result(MAXDIST2,1) - Result(MINDIST1,1); % Maximum amount that DIST1 can be shifted to ensure it is SOSD
tbase = Result(MINDIST2,1) - Result(MAXDIST1,1); % Maximum amount that DIST1 can be shifted back to ensure DIST2 is SOSD
% Initialize Golden Section Search upper and lower starting points
upper = 0;
lower = 0;
if tcomp > tbase
upper = tcomp;
lower = tbase;
elseif tcomp < tbase
upper = tbase;
lower = tcomp;
else
upper = tbase + 10;
lower = tbase - 10;
end
% Initialize Flags to Test for Convergence
bail = 0;
converge = 0;
while converge ~= 1 && bail < bailthresh % Find the minimum proportion that makes DIST1 SOSD DIST2
middle = (lower + upper) / 2;
if SOSDIntegralTestv3(DIST1 + middle, DIST2) == 1
upper = middle;
else
lower = middle;
end
if lower > upper
error('lower > upper!')
end
if (upper - lower) <= thresh % Convergence acheived when upper and lower are within thresh tolerance
converge = 1;
end
bail = bail + 1;
end
if converge == 1
Result(DIST1MINPROP,1) = upper;
if Result(MEANDIST2,1) > 0
Result(RELCOMPSOSDBASE,1) = upper / Result(MEANDIST2,1);
end
else % Golden Section Search Failed to converge because Threshold not met
Result(DIST1MINPROP,1) = -777777;
Result(RELCOMPSOSDBASE,1) = -777777;
end
% Initialize Golden Section Search upper and lower starting points
if tcomp > tbase
upper = tcomp;
lower = tbase;
elseif tcomp < tbase
upper = tbase;
lower = tcomp;
else
upper = tbase + 10;
lower = tbase - 10;
end
% Initialize Flags to Test for Convergence
bail = 0;
converge = 0;
while converge ~= 1 && bail < bailthresh % Find the minimum proportion that makes DIST2 SOSD DIST1
middle = (lower + upper) / 2;
if SOSDIntegralTestv3(DIST2, DIST1 + middle) == 1
lower = middle;
else
upper = middle;
end
if lower > upper
error('lower > upper!')
end
if (upper - lower) <= thresh % Convergence acheived when upper and lower are within thresh tolerance
converge = 1;
end
bail = bail + 1;
end
if converge == 1
Result(DIST2MINPROP,1) = lower;
if Result(MEANDIST2,1) > 0
Result(RELBASESOSDCOMP,1) = lower / Result(MEANDIST2,1);
end
else
Result(DIST2MINPROP,1) = -777777;
Result(RELBASESOSDCOMP,1) = -777777;
end
if Result(PRCROPFAILDIST1,1) >= 0 && Result(PRCROPFAILDIST2,1) >= 0
Result(CFDIFF,1) = Result(PRCROPFAILDIST1,1) - Result(PRCROPFAILDIST2,1);
else
Result(CFDIFF,1) = -999999;
end
% Categorize Risk: DIST1 More Risky (-1)/less Risky (1)/ Indeterminant (0) compared to DIST2
if Result(DIST1MINPROP,1) <= -777777 || Result(DIST2MINPROP,1) <= -777777
Result(DELTARISK,1) = -999999;
elseif Result(DIST1MINPROP,1) > 0 && Result(DIST2MINPROP,1) > 0
Result(DELTARISK,1) = -1;
elseif Result(DIST1MINPROP,1) < 0 && Result(DIST2MINPROP,1) < 0
Result(DELTARISK,1) = 1;
else
Result(DELTARISK,1) = 0;
end
end
2.3.0.4 CreateWTPBoundsbyCell
% INPUT Matrices: ndataseq01060216 and ndataseq10060216 from crop
% simulation models with Columns
% 1 = Unique Cell ID
% 2 = Replication
% 3 = Area
% 4 = Base
% 5 = New
%
%
% OUTPUT Matrix: RA with columns
% 1 = cell30m
% 2 = Comparison ID
% 3 = Base ID
% 4 = Mean Yield for Comp
% 5 = Standard Deviation of Yield for Comp
% 6 = CV of Yield for Comp
% 7 = Maximum Yield for Comp
% 8 = Minimum Yield for Comp
% 9 = Probability of Crop Failure for Comp
% 10 = Min Proportion for Comp to SOSD Base
% 11 = Mean Yield for Base
% 12 = Standard Deviation of Yield for Base
% 13 = CV of Yield for Base
% 14 = Maximum Yield for Base
% 15 = Minimum Yield for Base
% 16 = Probability of Crop Failure for Base
% 17 = Min Proportion for Base to SOSD Comp
% 18 = Difference in mean Comp - Base
% 19 = Difference in standard deviation Comp - Base
% 20 = Difference in CV Comp - Base
% 21 = Difference in Prob of Crop Failure Comp - Base
% 22 = Min Proportion for Comp to SOSD Base divided by average base yield
% 23 = Min Proportion for Base to SOSD Comp divided by average base yield
% 24 = Comp More Risky (-1)/less Risky (1)/ Indeterminant (0) compared to Base
% 25 = Area
% Compare Q to G
RA_yearseq_01 = ConstRiskThresh(ndataseq01060216, 4, 5);
LEN = length(RA_yearseq_01);
apones = ones(LEN,1);
clear apones LEN;
RA = [RA_yearseq_01];
RA=RA(1,:)
We can see here now that G stochastically dominates Q.
5 Bihar Model Application
The example above has shown for a case where there is one location in which there are multiple years of data. In the case of when one has run a crop growth model and has ncdf4 outputs, then one can use the same approach to estimate the model for the whole state.
Code
# Subset for each scenario and for Bihar and then send cell in row and var_year in columns# load some librarieslibrary(rasterVis)library(raster)library(rgdal)library(colorspace)library(RColorBrewer)#library(tidyverse)library(ggplot2)library(reshape2)library(ggtext)library(plyr)# define key files#load files to be used laterfile1 <-"data/simulations/s001.nc4"file2 <-"data/simulations/s002.nc4"file_aoi <-"data/aoi_mask/aoi_mask.tif"india_states <-"data/gadm40_IND_shp/gadm40_IND_1.shp"# loading key dataaoi <-raster(file_aoi)bihar <-readOGR(india_states)
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\MMKONDIWA\OneDrive - CIMMYT\Documents\GitHub\WTP_Bounds_SOSD_Risk_Model\data\gadm40_IND_shp\gadm40_IND_1.shp", layer: "gadm40_IND_1"
with 36 features
It has 11 fields
Code
aoi <- bihar[bihar$NAME_1=="Bihar",]plot(aoi)
Code
# functions for loading data load_sims <-function(files,var,aoi,operator,var2) {# check if operators is used, is NA, just load one varif (is.na(operator)==TRUE) { l <-list()#cycle through file names to load all different scenariosfor (i in files) { j <-which(i==files)#load file and add to list l[[j]] <- raster::stack(i,varname= var)[[3:35]] }} else { # if operator is used, load to vars with operators l <-list()for (i in files) { j <-which(i==files)# important: need to pass operator with backticks, not as stringif (var=="rice_yield") { a <- raster::stack(i,varname= var)[[3:35]] *3600}else {a <- raster::stack(i,varname= var)[[3:35]] }if (var=="wheat_yield") { b <- raster::stack(i,varname= var2)[[3:35]] *3340}else {b <- raster::stack(i,varname= var2)[[3:35]] } l[[j]] <-operator(a,b) } } l <-lapply(l,function(x){raster::mask(x,aoi)}) l <-lapply(l,function(x){raster::crop(x,aoi)})return(l)}# start extraction for ricecreate_csv <-function(crop,var) {files <-c(file1,file2)raster_df <-load_sims(files,var,aoi,NA)list_df <-lapply(raster_df, getValues)list_df <-lapply(list_df, na.omit)lapply(list_df, nrow)list_df <-lapply(list_df, "colnames<-", paste0(var,"_",1983:2015))#Scenarios# Fixed/onset is the rice planting strategy# Medium/long is the variety# suppl is where we ran supplementary irrigation instead of full irrigation. Use with caution.scen <-c("baseline","fixed_long")for (i in1:length(scen)) {write.csv(list_df[[i]],paste0("",crop,"_",scen[i],".csv"))}}create_csv("rice","rice_yield")create_csv("wheat","wheat_yield")
5.1 Clean for Octave
Code
# packages ---------------------------------------------------------------------library(dplyr)library(data.table)library(tidyr)library(rio)# RICE--------------------------------------------------------------------------## RICE-- Fixed long ---------------------------------------rice_fixed_long=read.csv("rice_fixed_long.csv")rice_fixed_long <-na_if(rice_fixed_long, '-99') #Replace all -99 with NArice_fixed_long <-na_if(rice_fixed_long, '-999999')rice_fixed_long=rename(rice_fixed_long,ID=X)rice_fixed_long$rice_area=rep(1,3386)rice_fixed_long_long <-melt(setDT(rice_fixed_long), id.vars =c("ID","rice_area"),value.name=c("fixed_long_rice_yield") ,variable.name ="year")rice_fixed_long_long$year=as.numeric(gsub('[^0-9]', '',rice_fixed_long_long$year))## RICE-- Baseline ---------------------------rice_baseline=read.csv("rice_baseline.csv")rice_baseline <-na_if(rice_baseline, '-99') #Replace all -99 with NArice_baseline <-na_if(rice_baseline, '-999999')rice_baseline=rename(rice_baseline,ID=X)rice_baseline_long <-melt(setDT(rice_baseline), id.vars =c("ID"),value.name=c("baseline_rice_yield") ,variable.name ="year")rice_baseline_long$year=as.numeric(gsub('[^0-9]', '',rice_baseline_long$year))rice_baseline_fixedlong=merge(rice_fixed_long_long,rice_baseline_long,by=c("ID","year"))rice_baseline_fixedlong$year=rice_baseline_fixedlong$year-1982rice_baseline_fixedlong=as.data.frame(rice_baseline_fixedlong)rice_baseline_fixedlong$baseline_rice_yield[is.na(rice_baseline_fixedlong$baseline_rice_yield)]=0rice_baseline_fixedlong$fixed_long_rice_yield[is.na(rice_baseline_fixedlong$fixed_long_rice_yield)]=0rice_baseline_fixedlong=rice_baseline_fixedlong[,c(1,2,3,5,4)]export(rice_baseline_fixedlong,colNames=F,"rice_baseline_fixedlong_s.xlsx")library(rhdf5)rice_baseline_fixedlong=as.matrix(rice_baseline_fixedlong)colnames(rice_baseline_fixedlong) <-c(1,2,3,4,5)h5save(rice_baseline_fixedlong,file="rice_baseline_fixedlong.hdf",createnewfile =TRUE)rice_baseline_fixedlong=as.data.frame(rice_baseline_fixedlong)# Remove all zeros and nasrice_baseline_fixedlong_nonzero=subset(rice_baseline_fixedlong, rice_baseline_fixedlong$fixed_long_rice_yield!=0)rice_baseline_fixedlong_nonzero=subset(rice_baseline_fixedlong_nonzero, rice_baseline_fixedlong_nonzero$baseline_rice_yield!=0)export(rice_baseline_fixedlong_nonzero,colNames=F,"rice_baseline_fixedlong_nonzero.xlsx")library(rhdf5)rice_baseline_fixedlong_nonzero=as.matrix(rice_baseline_fixedlong_nonzero)colnames(rice_baseline_fixedlong_nonzero) <-c(1,2,3,4,5)h5save(rice_baseline_fixedlong_nonzero,file="rice_baseline_fixedlong_nonzero.hdf",createnewfile =TRUE)# WHEAT ---------------------------------------------------------------------------------------------------## wheat-- fixedlong ---------------------------------------wheat_fixedlong=read.csv("wheat_fixed_long.csv")wheat_fixedlong <-na_if(wheat_fixedlong, '-99') #Replace all -99 with NAwheat_fixedlong <-na_if(wheat_fixedlong, '-999999') #Replace all -99 with NAwheat_fixedlong=rename(wheat_fixedlong,ID=X)wheat_fixedlong$wheat_area=rep(1,3386)wheat_fixedlong_long <-melt(setDT(wheat_fixedlong), id.vars =c("ID","wheat_area"),value.name=c("fixedlong_wheat_yield") ,variable.name ="year")wheat_fixedlong_long$year=as.numeric(gsub('[^0-9]', '',wheat_fixedlong_long$year))## wheat-- Baseline ---------------------------wheat_baseline=read.csv("wheat_baseline.csv")wheat_baseline <-na_if(wheat_baseline, '-99') #Replace all -99 with NAwheat_baseline <-na_if(wheat_baseline, '-999999') #Replace all -99 with NAwheat_baseline=rename(wheat_baseline,ID=X)wheat_baseline_long <-melt(setDT(wheat_baseline), id.vars =c("ID"),value.name=c("baseline_wheat_yield") ,variable.name ="year")wheat_baseline_long$year=as.numeric(gsub('[^0-9]', '',wheat_baseline_long$year))wheat_baseline_long=merge(wheat_fixedlong_long,wheat_baseline_long,by=c("ID","year"))wheat_baseline_long$year=wheat_baseline_long$year-1982wheat_baseline_long=as.data.frame(wheat_baseline_long)wheat_baseline_long$baseline_wheat_yield[is.na(wheat_baseline_long$baseline_wheat_yield)]=0wheat_baseline_long$baseline_wheat_yield[is.na(wheat_baseline_long$baseline_wheat_yield)]=0wheat_baseline_long=wheat_baseline_long[,c(1,2,3,5,4)]export(wheat_baseline_long,colNames=F,"wheat_baseline_long_s.xlsx")library(rhdf5)wheat_baseline_long=as.matrix(wheat_baseline_long)colnames(wheat_baseline_long) <-c(1,2,3,4,5)h5save(wheat_baseline_long,file="wheat_baseline_long.hdf",createnewfile =TRUE)# Remove all zeros and naswheat_baseline_long=as.data.frame(wheat_baseline_long)wheat_baseline_long_nonzero=subset(wheat_baseline_long, wheat_baseline_long$baseline_wheat_yield!=0)wheat_baseline_long_nonzero=subset(wheat_baseline_long_nonzero, wheat_baseline_long_nonzero$baseline_wheat_yield!=0)export(wheat_baseline_long_nonzero,colNames=F,"wheat_baseline_long_nonzero.xlsx")library(rhdf5)wheat_baseline_long_nonzero=as.matrix(wheat_baseline_long_nonzero)colnames(wheat_baseline_long_nonzero) <-c(1,2,3,4,5)h5save(wheat_baseline_long_nonzero,file="wheat_baseline_long_nonzero.hdf",createnewfile =TRUE)h5closeAll()
